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Monolayer graphene at neutrality in the quantum Hall regime has many competing ground states 
with various types of ordering. The outcome of this competition is modified by the presence of 
the sample boundaries. In this paper we use a Hartree-Fock treatment of the electronic correla¬ 
tions allowing for space-dependent ordering. The armchair edge influence is modeled by a simple 
perturbative effective magnetic field in valley space. We find that all phases found in the bulk of 
the sample, ferromagnetic, canted antiferromagnetic, charge-density wave and Kekule distortion are 
smoothly connected to a Kekule-distorted edge. The single-particle excitations are computed taking 
into account the spatial variation of the order parameters. An eventual metal-insulator transition 
as a function of the Zeeman energy is not simply related to the type of bulk order. 


I. INTRODUCTION 

When subjected to a strong perpendicular magnetic field, the electrons confined to the two-dimensional (2D) carbon 
lattice of graphene form a unique quantum Hall (QH) system. Notably graphene at neutrality is an example of QH 
ferromagnetism with many competing ground states!"**. j n g U ch systems we expect generically complex physics at 
the edge. Early work on graphene edge stateiP^ has shown that when taking into account the electron spin degree 
of freedom, the edge states should be of helical nature, i.e., exhibiting counterpropagating modes carrying opposite 
spin polarizations. Graphene was thus proposed to be a model candidate for a quantum spin Hall (QSH) system. 
Furthermore, graphene has unique features so that it may be an ideal probe material to study QH edge physics 
experimentally : in semiconductor-based 2D electron gas systems the confining potential is soft at the scale of the 
magnetic length so that there may be edge reconstructiorP^utEl, spoilin g the ideal theoretical description. Therefore, 
the understanding of edge phenomena in these 2D systems is difficult® 16 ' and still controversial. In contrast, the 
boundaries of the graphene lattice are directly defined at the atomic level. Therefore, they naturally represent 
atomically sh arp Q H edges. This should allow observation of QH edge state physics without complications from edge 
reconstruction 16 11 . The fabrication, design, and control of graphene edge structures with atomic level precision is a 
field of ongoing research : see, e.g., Refs. [T8H201 Among the theoretical approaches, also a mean-field treatment of 
a Hubbard-type model has been applied to the edge physic^. The QSH nature of the graphene edge has been the 
subject of recent experimental investigation:^. At Landau level filling factor v = 0, upon tilting the applied magnetic 
field with respect to the graphene sheet, there is a metal-insulator transition for some critical angle. This suggests a 
change of the bulk state as a function of tilting. Indeed, extensive theoretical studies of the v = 0 ground state (GS) 
structure of bulk graphen d 6 F I 2 :j l 24 l have shown the existence of various competing phases with distinct symmetry¬ 
breaking properties. While the graphene GS at neutrality is a highly degenerate SU(4) ferromagnetic multiplet, small 
symmetry-breaking terms due to short-distance physics lift this degeneracy and the system may form various different 
phases. Among these phases are e.g. the ferromagnet (F) statd^or the antiferromagnet (AF) state, where the latter 
may form a canted antiferromagneti c (CAF) state as has been pointed out by Kharitonov in Ref. [7J Further possible 
phases include a charge-density waW- 5 - 5 Es| 26 | (CDW) or a Kekule distorted stat d 27 * 2 ^ (KD). Transitions between these 
states may be induced by varying the Zeeman energy which is done by tilting the field. 

In this paper we study the edge structure of the v = 0 QH state in the presence of SU(4)-symmetry-breaking 
interactions. We use a simple model of the edge potential in the basis of the n = 0 LL orbitals appropriate to an 
armchair termination of the graphene lattice and treat interactions by a Hartree-Fock (HF) approximation. Our 
variational ansatz is orbital-dependent so it captures the spatial variations of the ordering from the bulk to the edge 
(an effect which is absent in previous HF studied^. We find that there is always a crossover towards a Kekule 
distorted region close to the edge. There is appearance also of spin/valley nontrivial entanglement which is limited 
to the transition region towards Kekule order and does not take place either in the bulk or close to the edge. We 
propose a quantitative measurement of the entanglement by computing the concurrence as a function of edge distance. 
Within HF we also compute the single-particle properties of the particle-hole excitation spectrum. We discuss how this 
spectrum varies with the edge distance and also how it is influenced by the nature of the bulk order. Our main finding 
is that there should be always a metal-insulator transition as a function of the Zeeman field whatever the nature of the 
ordered phase. So strictly speaking the experimental observations of Ref. [22] do not imply that their graphene bulk 
is a CAF state. It should be noted that our HF treatment has some shortcomings. Notably the Coulomb exchange 
interaction in the n = 0 LL leads to a coupling between charge and spin/valley degrees of freedom. So in general 
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charge motion is done through spin/valley textures as proposed in Refs. fTTTThfTITTl As in previous HF workPH we will 
not try to model these effects. While the edge effects will overcome exchange energy close enough to boundary, they 
may change the nature of excitations right in the transition region. More work is needed to understand this point. 

The paper is structured as follows : In Sec. [TTJ we define the theoretical framework describing the v = 0 QH 
state of graphene in the presence of a boundary and the corres pond ing m odel Hamiltonian. We introduce the 
parametrization of the Hartree-Fock (HF) GS wave function in Sec. Ill In Sec. IV, we present results for the GS and 


its properties obtained from minimizing the HF energy functional. We describe the evolution of the different possible 
bulk phases when moving towards to the edge. We find that the presence of a boundary gives rise to novel spin/isospin 
configurations which do not exist in the bulk. Furthermore, we show the existence of an intermediate region between 
the bulk and the edge with nontrivial entanglement of spin and valley isospin. In Sec. [VJ we extend our HF treatment 
to compute the spatial evolution of the single particle (SP) energy levels and corresponding eigenstates from the bulk 
to the boundary. As pointed out by previous worlPS, the SP spectra can either exhibit nonzero gaps or support 
gapless edge states, depending on the system parameters determining the bulk phase. We describe how the spatial 
variation of the spin/isospin texture influences the shape of the SP energy levels. This leads to an understanding of 
the edge gap, the number of conducting channels, as well as possible conclusions for the bulk symmetry properties 
drawn from the conductance behavior of the edges. In a final part, we compute the spin and isospin properties of the 
corresponding SP eigenstates to directly prove that the edge states of the v = 0 QH state in graphene indeed exhibit 
the helical properties of a QSH state. Sec. |VI| finally discusses our results in relation to experimental findings and 
contains our conclusions. 


II. MODEL OF THE GRAPHENE EDGE 


We first recall basic facts about the electronic structure of graphen d' 34 ! 35 ! The hexagonal structure admits two 
triangular Bravais sublattices A and B that form the basis for a tight-binding Hamiltonian. In the Brillouin zone 
there are two special degeneracy points : at these Dirac points K and K ', the valence and the conductance band form 
Dirac cones and touch at the Fermi level for neutral graphene. In the vicinity of the Dirac points, for each orientation 
of the spin a =t,|, the electronic wave functions 'Fa and on the two sublattices can be written as : 


'pA,<r(r) = e lK ' r V’A, CT + e ,K ' T ip ' A)(7 , 

'f'B,o-(r) = e iK ' r ipB,a + e“* K ' T ^' B ,a- 

Assembling the envelope functions in a four spinor notation we write for the electronic state : 
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( 2 ) 


where the subindex Hk,k' ® Ha,b indicates that the state H/o- lives in the Hilbert space formed as the direct product 
between Dirac valley space Hk,k' and the A, B sublattice space Ba,b■ An applied magnetic field leads to the 
formation of Landau levels (LL) with energies : 

E n =^~ \/2 \ti\ sgn(n), (3) 


where n £ Z is the LL index, £b = y is the magnetic length for perpendicular magnetic field strength B ±, and 
Vp denotes the Fermi velocity. The corresponding eigenstates can be written as : 
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The filling factor v of the Landau levels is defined by : 
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FIG. 1: Shape of the kinetic energy edge potential Ekin in T-Lkin of Eq. |8|. The kinetic energy rises from Ekin = 0 in the bulk 
(equivalent to the case of infinitely extended, translationally invariant graphene in the n = 0 LL) to the energy of the n = 1 LL 
exactly at the edge. The curve was obtained from solving the problem of free electrons on a graphene lattice in the presence 
of a magnetic field in the presence of boundary, i.e., by applying appropriate boundary conditions for an armchair edge (see 
inset). Zigzag edges have similar nonzero energy states but also additional zero-energy modes that are beoynd the scope of our 
work. 


where n e denotes the electronic density. The configuration of neutral graphene, i.e., v = 0 is peculiar. Indeed this 
particle-hole symmetric situation corresponds to the case in which all the LLs with n < 0 are filled and all the LLs 
with n > 0 are empty, but the n = 0 LL is exactly half filled with two electrons per orbital. The form of the n = 0 
wave function as given in Eq. Q implies that n = 0 LL electrons reside on one of the sub-lattices only. In the 
following, we consider the case of neutral graphene and therefore study the properties of electrons in the n = 0 LL. 
We simplify the notation by collecting only the non-zero entries of the n = 0 spinor in Eq. @ as : 


T n = 



( 6 ) 


identifying the valley and the sublattice indices in a common valley-isospin rasT = + = K<->A and r = — = 
K' o B. In the four-dimensional Hilbert space H = H sp i n ® H va u ey we use the indices /x, u to label the four possible 
configurations of spin and isospin in this space: /i, v £ {t +> t — > i +, I —}. Due to the fourfold degeneracy in spin 
space (a =|, |) and in valley space (r = +, —), integer QH effect in graphene is expected at values of v that change 
in steps of four. 

The total Hamiltonian we use is given by : 


R = Rkin + R-Z + Rcoul + RAnisoi (?) 

where we have a space dependent kinetic energy induced by the presence of the boundary : 

Rkin = - ^2 E kin(^i)rl, ( 8 ) 

i 

where the index i labels the positions of the electron orbits. The hexagonal graphene lattice can be terminated in 
many different ways, yielding several possible edge structures. Every different atomic configuration leads to different 
boundary conditions for the wavefunction 36 . Two extreme cases are the so-called zigzag and armchair edges 30 . A 
finite piece of graphene terminated by a zigzag edge and an armchair ed ge is show n in the inset of Fig. [l] The kinetic 
energy and the corresponding eigenstates can be obtained analyticall ^ 10 l 31 l 3 lH8|. This is equivalent to turning the 
level index into a space dependent quantity n(R) where R relates to the distance to the edge r as r = R£b/V 2. 
In Fig. [T] we show the spatial shape of the kinetic energy Ekin obtained by this procedure, as we will use it in the 
subsequent calculations. We write the kinetic energy as a space-dependent potential proportional to t x (a comparable 
treatment can be found in Refs. I3SI53) . This corresponds to a perturbative treatment as it assumes an expansion of 
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FIG. 2: Bulk GS phase diagram as a function of the coupling energies u± and u z of the v = 0 QH state for a Hamiltonian as in 
Eq. Q. Different colors distinguish between the different possible GS phases. The white, dotted lines indicate the parameters 
we choose throughout this paper to perform cuts through the phase diagram in order to explore the behavior of all possible GS 
phases when starting from the bulk and moving towards an edge of the graphene lattice. The phase diagram for the different 
GS phases in the bulk of graphene was first presented in Ref. [7j 


the perturbed edge states in terms of the unperturbed bulk basis states. It restricts our description to the case of 
” armchair-like” boundaries : one can always infer the number of branches in the SP edge spectrum as being equal to 
the number of degenerate SP levels in the bulk and hence apply a perturbative expansion as implied by Eq. ([8]). A 
derivation of such a Hamiltonian describing the kinetic potential of a graphene edge using arguments of perturbation 
theor y can be found in Ref. 1291 The edges terminated by a zigzag boundary condition support dispersionless surface 
states 30 31 that break the simple bulk/edge correspondence. They are beyond our simple treatment. The form of the 
kinetic energy in Eq. (8]) is valid only in the regime E^ in <C hui c , i.e., spatially not too close to the edge. As can 
be seen from Fig. [l] this condition is very well met if we restrict the subsequent discussion to the regime R > 3. 
Hence the restriction R > 3 corresponds to a minimal distance r m i n ss 2.121b, which at realistic experimental values 
corresponds to r m j„ « 120ao, where ao denotes the lattice constant of graphene. The Zeeman energy can be written 
as : 


Hz 



(9) 


In the case of graphene the spacing between kinetic LL energy levels, A Ek in , easily exceeds the Zeeman energy by 
two orders of magnitude. The Coulomb interaction is given by : 


R-Coul 




( 10 ) 


where e is an effective dielectric constant which depends upon the substrat^H. It has full SU(4) symmetry. At 
neutrality we have an example of quantum Hall ferromagnetisirfD with this large symmetry : the ground state is thus 
highly degenerate and forms an irreducible representation of SU(4). However this symmetry is only approximate. In 
fact it is weakly broken by lattice-scale effects that include short-range Coulomb interactions and electron-phonon 
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couplings. It is difficult to obtain precise estimates of these effects but their symmetry-breaking properties can be 
encoded in the following Hamiltonian : 

T-LAniso = [fi , *' r * r * + 9yT l yT. J y + g z r l z T J z ] S 2 {Ti - rj), (11) 

with <5 denoting the Dirac delta function. This Hamiltonian 'H.Aniso has been proposed by Aleiner et a/P^and its effects 
have been analyzed at the mean-field level by Kharitonov 7 . Its symmetry properties and phase diagram have been 
studied by exact diagonalization 41 . It is parametrized by the coupling constants g x ,y,z whose values are not known 
with precision. It is likely that the ratio of the energy scales between Coulomb interaction and these anisotropies is 
of the order of 10 2 . It is thus best to explore the complete phase diagram taking these parameters as unkn owns . For 
monolayer and bilayer graphene at neutrality, there is a rich phase diagram as a function of these coupling ;:l 7 f 12 ( The 
fractional quantum Hall states are also sensitive to these effect 

We now perform a HF study of this Hamiltonian by including the edge potential. We note that in this approach 
we neglect all possible spatial dependence of the coupling constants, which is justified as long as we analyze a spatial 
domain not too close to the edge. 


III. HARTREE-FOCK TREATMENT 


The neutral v = 0 state corresponds to the half-filled case where two of the four available states per orbital are 
occupied. We look for the ground state within the family of Slater determinants of the form : 


|G) = II ( SW cl(p)ct(p) I |0>, 


( 12 ) 


where p denotes the Landau-gauge momentum component along the edge which labels the electron orbitals. The 
vacuum |0) consists of the completely occupied set of states for all n < 0 and completely empty states for all n > 0. 
In Eq. (121 g is a 4 by 4 antisymmetric matrix, i.e., g pl/ = —g vli , in order to describe a valid fermionic state and 
Tr[f/f/'] = 2 to ensure normalization of the two-particle state. We minimize the energy of the Slater determinant by 
varying g. To capture the effect of the edge potential we take the g matrix to be momentum dependent, i.e., g = g{p) 


in Eq. (12). Due to the duality between the longitudinal momentum p and the transverse coordinate r p = p£ 2 B this 
is equivalent to a space dependent description of the problem. The most general antisymmetric matrix g has 12 real 
parameters. By exploiting the symmetry properties of the state |G) and the Hamiltonian, one can reduce the number 
of free parameters. We use the same strategy as in Ref. where an equivalent problem was studied in the context 
of electronic bilayer systems. It is convenient to rewrite the problem in terms of the following simple expectation 
values : 


S a = ^(G\cH P )a a c(p)\G) = ^Tr[a a gg^], 

T a = ^ (G|c f (p)T a c(p) | G) = ^Tr [r aff5 t ], 

Rap = ^(G\J (p)a a rpc(p)\G) = ^Tr[cr Q Tagp t ]. 


(13a) 

(13b) 

(13c) 


The expressions of Eqs. (13a) and (13b) yield the components of the total spin S a and isospin T a per orbital p. The 


energy contribution from the symmetry breaking interaction Hamiltonian of Eq. (Ill is now given by : 

{G\n A ni SO \G) = ^ ^2 u <* ( T - TrlTaSffVaffS 1 ]) , ( 14 ) 

a 

where the anisotropy energies u a are given by u a = gff t- Isotropy of the interaction potential in the x-y-plane 


implies that u x = u y =: u±. Using Eqs. (13) and (14), we obtain the following expression for the functional of the 
total energy E tot = {G\H\G) : 


Etot = -2 E kin T x - 2 E z S z + ^ u a (t 2 - ^ Rj a - S 2 ). 


(15) 


The 12 — 2 = 10 free parameters of the problem (dropping the overall phase and normalization constant) are now 
encoded in the 6 components of the total spin S and the total isospin T, together with 4 out of 9 components of R a p 
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which can be chosen independently. The invariance of E to t in Eq. (15) under rotations of S in spin space and rotations 
of T around the 2 -axis in isospin space allows us to choose S y = S x = T y = 0 with no loss of generality, yielding 
seven variables to be determined. The dimension of parameter space can be further reduced by careful consideration 
of all the symmetries of the problem. As demonstrated by Ezawa et oiP^O in a situation of an equivalent symmetry 
class, reduction is possible to a total number of three free parameters. For the present system, this leads us to a 
minimization problem for the total energy E tot with respect to the variational parameters — 1 < a < 1,-1 < /3 < 1, 
and X S IR, which are related to observables of Eq. (131 by : 


S x = , 1 = \/1 — a 2 , T x = — a\/l — /3 2 , T z = u X = ap, 


X 2 


\/r+x 2 


yr+x 2 


(16) 


and 


E n>2 _ ±z_ ST' f>2 2o2 n2 ± x 

^.2 ’ Z_^ ° > Z-* Wz ^,2 ’ 


(17) 


where the index i runs over the spatial components {x,y, z}. The density matrix p g = gg 1 is connected to these 
quantities as (summation convention implied) : 


Pg — 2 “ ^ 2 ff 7 T R 7v T u,Ty Rij'j • 


(18) 


IV. GROUND STATE PROPERTIES 


A. Evolution of the Spin-Isospin Texture close to the Edge 


The bulk GS of graphene within the model we use has several different phases depending on the anisotropy energies 
u_i and u z compared to the Zeeman energy E ^ 41 [ The anisotropies u±, u z and the Zeeman term Ez select some subset 
of the manifold of SU(4) ferromagnetic ground states. These phases have distinct spin and isospin configurations, 
i.e., by different spin textures. The mean-field GS diagram is shown in Fig. [2] It is correct beyond mean-field as 
shown by exact diagonalization technique:!^. Notably all these bulk phases do not involve spin/valley entanglement. 
We now generalize the description of these phases by treating the effect of a boundary of the lattice. We proceed as 
follows : we minimize the energy E to t(ct, P,x) °f Eq. ( fl5| , including a space dependent edge potential of the shape 
shown in Fig. [l] by varying the parameters a,/3,x- Then from the knowledge of the parameters a(R), /3{R), x(-R)> 
we can compute the values of the observables S Z (R),T X (R),T Z (R) of the GS |G) via Eq. (161. It is also possible to 
construct the entire density matrix p g characterizing the GS via Eq. (18). 

In order to obtain a full picture capturing the edge behavior of all possible bulk phases shown in Fig. [2j our choice 
of system parameters is guided by the following idea : for fixed Zeeman energy Ez, we vary the coupling energies 
u_ l and u z because this can be realized experimentally by tilting the magnetic field. We choose three values of the 
perpendicular coupling energy u z : u z = 5Ez, u z = 2 Ez and u z = —2 Ez, and we vary the perpendicular coupling in 
the range —3 Ez <u _l < 3 Ez- This leads to horizontal cuts through the v — 0 GS phase diagram in the u_l, u 2 -plane, 
shown by white dotted lines in the phase diagram of Fig. [T] For u z = 5Ez and u z = 2 Ez by varying u± we meet the 
KD, CAF and F phases : 

KD . CAF . F 

u± = — oo -1 - 1 --- > u± 


u± > - \{u z + y/2E% 


+ u 2 z ) 


u± > -if 


For u z = —2 Ez and varying again u± we find the KD, CDW and F phases : 
KD . CDW 


U _L = —OO 


U _L > U z 


u_l > -(Ez + u z ) 


-> U _L 


The corresponding bulk phase transitions are indicated by white arrows in the phase diagram in Fig. [I] 

We first investigate the influence of the edge potential on the spin and isospin observables S and T. More precisely, 
we discuss the spatial evolution of the components S z ( R ), T x (R ), T z (R) for different choices of the anisotropy energies 
u± and u z compared to the Zeeman energy Ez- Figure [3] shows the results for u z = 5Ez , corresponding to a cut 
through the upper half plane (we plot obsevables with colored lines). The four different panels depict the situation 
for different values of u± : for u± = —Ez, the bulk system is in the CAF phase with canting angle cos 9 = S z = 1/2 
(upper left panel), whereas for u± = —0.2 Ez, u± = 0.5 Ez, and u± = 1.5 Ez the bulk system establishes a F phase in 
which the spins are fully polarized. The colored curves shown in Fig. [4] correspond to the obervables for a cut through 
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FIG. 3: Evolution of the spin and isospin as well as the concurrence C as functions of R = \f2r/lB , with r the distance from 
the edge. We fix u z = 5Ez and vary u±. (Colored lines: S z , solid; T z , dashed; T x , dotted. Black, thin, solid line: Concurrence 
C). There is a transition between the bulk phase on the right-hand side [CAF for u± = —Ez in the upper left panel (green 
colors) and F for all other choices of u± shown (blue colors)] to a KD phase at the edge. In an intermediate regime spin and 
isospin, are canted with 0 < S z < 1 and 0 < T x < 0 at the same time and nonzero values of the concurrence. With growing 
Mj_, this domain wall grows narrower in space and moves closer to the boundary. 



FIG. 4: Same as Fig. [3] for transverse coupling energy u z = —2 Ez and different perpendicular coupling energies u_l, such that 
the phases of the lower half plane of the GS phase diagram are established in the bulk. Left panel: u±_ = 0.6, leading to a 
CDW in the bulk (red colors). Right panel: u± = 1.2 Ez, at which the bulk is in a F phase (blue colors). 


the lower half plane of the phase diagram at u z = —2 Ez- Here, the perpendicular couplings are chosen so that the 
left panel at u± = —0 .6Ez corresponds to a CDW phase whereas the right panel at u± = 1.2 Ez again corresponds to 
a F bulk phase, as predicted by the GS phase diagram of Fig. [2j Curves for values of the anisotropy energies favoring 
a KD phase in the bulk are not shown since in this case the system does not undergo any transition whatsoever but 
remains in the bulk KD phase all the way to the edge. 

In general, one can distinguish three different regimes for the behavior of the observables as a function of the 
distance r = ^R to the edge. For sufficiently large values of R, i.e., deep enough in the bulk, we recover the results 

of mean-field theory 7 . Close enough to the edge, the system is driven into a KD phase with N x = 1 and S z = N z = 0, 
independently of the bulk phase it adopts. This behavior is due to the edge potential in the kinetic energy Hamiltonian 
Rkin hi Eq. ([8]) : this term is proportional to t x , so it acts as a Zeeman effect in isospin space, polarizing the isospin 
along the ^-direction as soon as Ekin(R) is large enough. This behavior is also consistent with previous workJ^EU. 
In an intermediate regime we find a finite interval in space in which S z ^ 1, T x ^ 1 and T z ^ 0, N x ^ 0, i.e., the 
spin and the isospin are canted simultaneously with respect to their bulk values. There is thus a domain wall at 
a small finite distance from the edge. For the CAF configuration, this domain wall connects smoothly to the bulk 
configuration. For a system in a F phase in the bulk however, the change in spin and isospin is abrupt and the domain 
wall is narrower with increasing u±. Hence, for larger values of zij_, the F phase of the bulk proves to be more resistant 
against the increasing influence of the edge. 

From our analysis and the results shown in Figs. [3] and [4] we therefore draw the following conclusions: The phases 

















































FIG. 5: Maximum concurrence C ma x in the different regimes from the bulk to the edge. For three different values of the 
coupling energy u z , we vary u±. Different symbols represent different phases in the bulk at the respective value of ui: F (blue 
circles), CAF (green squares), CDW (red diamonds), and KD (gray triangles). Empty, shaded, or filled symbols distinguish 
between the cuts at different u z . The dashed lines connect the data points as a guide to the eye. Vertical, gray lines mark the 
values of {u z ,u±}, at which, in the bulk, transitions between the respective phases occur. We observe that non-zero values of 
C happen only if the system in the bulk is in a CAF or in a F phase. In this case, the CAF/F transition is smooth. If the 
bulk is CDW or KD,the concurrence remains equal to zero all the way from the bulk to the edge. This leads to a discontinuous 
jump in the curve of C max at the point of the CDW/F transition. 


in the bulk of a finite sample of graphene do not remain unaffected close enough to the edges. Indeed the effective 
edge potential causes the bulk state to undergo a transition in which the polarizations of spin and isospin change 
simultaneously. Sufficiently close to the edge, the GS is driven into a KD phase independently of the nature of the 
bulk phase. 


B. Spin-Valley Entanglement of Edge States 


The parametrization of S,T, and R Qi g in terms of the three parameters a,/3, and x i n Eqs. (16) and allows 
for complete reconstruction of the density matrix p g via Eq. (18) using the values of the parameters obtained by 
minimizing E to t from Eq. (15). Therefore, we have a full description of the GS and its spatial evolution from the bulk 


towards the edge. In this paragraph we investigate the entanglement between spin and isospin degrees of freedom in 
the system. For the infinite bulk case, product states of the form |s) ® |n), where |s) denotes the (single p article) 
spin state and | n) the (single particle) isospin state, have been used as an ansatz to minimize_the GS energ}k d 29 l 42 l. 
Existing studies of edge states using a variational trial wave function approach have suggested^, however, that for a 
non-zero edge potential, spin and isospin might not remain independent, separable observables, but become entangled. 
In order to quantify the amount of entanglement in the bipartite two-level system H = H spin (g H va u ey , we calculate 
the concurrence C according to the defmitiorp3: 


C = max(Ai — A 2 — A 3 — A 4 ,0), (19) 

where the A * are the eigenvalues of the matrix 

= \fPg( (J y ® <J y)Pg{ <J y ® °y)\fPgi (20) 

in decreasing order Af > Af +1 Vi. In Eq. ( |20| ), o y denotes the 2x2 Pauli matrix. The quantity C ranges from 0 to 1 
with C = 0 meaning no entanglement ana (J = 1 for maximally entangled states. 

In order to study the entanglement of all phases, we perform the same cuts through the phase diagram as in the 
previous paragraph, fixing the transverse coupling at u z = 5Ez and u z = 2 Ez to investigate the upper half plane 
and at u z = — 2 Ez to study the lower half plane and we vary the perpendicular coupling u± at these fixed values. 
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FIG. 6: SP energy levels e±± for the different phase regimes as they are obtained from the analytical formulas in Table [T] 
Upper left panel: CAF (solid, green lines), F (blue, dashed lines). Upper right panel: CDW for |A| > Ez (solid, red lines) and 
|A| < Ez (dashed, black lines). Lower left panel: KD for |A| < Ez with A > 0 (solid, orange lines) or A < 0 (dashed, black 
lines). Lower right panel: KD for |A| > Ez with A > 0 (solid, orange lines) or A < 0 (dashed, black lines). Depending oil the 
sign and the magnitude of A, the different cases for the SP spectra differ in the number of level crossings, even within one and 
the same phase. 


Examples of the spatial behavior of the concurrence C{R ) as a function of the distance from the edge is depicted by 
the black solid lines in Figs. [3] and [4] The curves reveal several characteristics of the behavior of the concurrence. 
It goes to zero deep enough in the bulk for all values of the anisotropies liniR-K*, C(R) = 0 V u z ,u±. Close enough 
to the edge, the concurrence is also equal to zero for all possible bulk phases, as can be seen in Figs. [3] and [I] for 
C(R ~ 3) = 0 V u z ,u±. In an intermediate regime for which the system is in a F or CAF phase in the bulk, we find 
that the concurrence develops a sharp peak. This peak appears precisely within the domain wall separating its bulk 
phase to the KD phase near the edge. The peak is sharper and higher with rising u± , as the domain wall becomes 
more and more narrow in space. Another behavior is observed when the bulk is CDW (left panel of Fig. [4]). Here, 
the concurrence remains zero independently of the distance from the edge : C(R) = 0 \/R. 

Our findings are summarized in Fig. [5j where we plot the maximum concurrence C max as a function of u±. The 
resulting curves characterize the behavior of the spin-valley entanglement of the edge states. Nonzero values of the 
concurrence are found only for anisotropies favoring a CAF phase (green squares) or a F phase (blue circles) in 
the bulk. In these regimes, the maximum concurrence C max is a monotonically rising function of the perpendicular 
coupling itj_, with no discontinuity at u± = —Ez/ 2, which would correspond to the CAF/F transition in the bulk. 
Discontinuous jumps appear at values of «j_ corresponding to the transitions KD/CAF or CDW/F. Combining the 
information from Figs. ii and Fig. [5] we draw the following conclusions : unlike the states in an infinite system, 
the GS in the presence of a boundary may exhibit nonzero spin-valley entanglement. As demonstrated in Figs. [3] and 
[4j the concurrence is exactly zero in all configurations where either the spin or the isospin is strictly zero. Nonzero 
values of the concurrence appear for configurations in which both spin and isospin are canted simultaneously. 

Compared to the bulk case, we find that the lattice boundary gives rise to novel ground state phases close to the 
edge, with simultaneous canting of spin and isospin, 0 < S < 1, 0 < T < 1, and nonzero spin-vall ey e ntanglement. 
These phases cannot be described using trial wave functions in the form of separable product states 7 ^. 


V. MEAN FIELD SPECTRUM OF EXCITED STATES 
A. Single Particle Energy Levels 

We now diagonalize the single-particle HF Hamiltonian. The spectrum of the excited states is of particular interest, 
since the conduction properties of real graphene samples are governed by the edge modes in the QH regime. Recent 
conductance experiments have shown 22 that upon tilting the applied magnetic field there is a transition from an 
insulating regime to a phase where presumably edge states carry a nonzero current. Tilting the magnetic field 
corresponds to varying the parameter u±/Ez in the system. The experimental observations therefore suggest that 
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the gap to excited states in the edge spectrum varies as a function of Ui_/Ez and closes, eventually, giving rise to a 


metal-insulator transition. We write the one-body HF Hamiltonian h 
Eq. Q. It consists of four terms : 


HF 


corresponding to the full Hamiltonian T~L of 


Kv (P) = ~Ekin(p)[Tx\»iy ^ Ez[^z]imu + A Ml/ + ' 4 A 




( 21 ) 


which we obtain via standard HF decoupling. For the mean-field potential from the Coulomb interaction Hamiltonian 
of Eq. (10 1 we find 


= -u 0 [gg ] ] h 


( 22 ) 


where uq describes the exchange term of the Coulomb interaction Hamiltonian of Eq. (10). This formula is valid 
provided we neglect the spatial dependence of g. It means that we do not capture the spin texture effects of the 
exchange. For completeness, in the following analytical calculations and expressions, the Coulomb contribution of 


Eq. (22) will be written explicitly. The mean-field potential due to the interactions breaking SU(4) symmetry is given 
by : 


a y ) u c 


{y a ]^^[gg^T a } - [T a gg^i 


(23) 


Within HF mean-field theory, diagonalizing li provides access to SP energies £», satisfying Ir |z) = £i\i), where 
the state labeled by i stands for the zth SP HF eigenstate. In the following, we assume the eigenvalues to be ordered 
£i < £2 < £3 < £ 4 . In earlier work 2 ^, this Hamiltonian has been studied with the assumption of constant order 
parameter to the edge. 

However, the explicit effective valley field due to the edge certainly invalidates this simple assumption. It should 
be noted also that in the intermediate regime between the bulk state and the edge regime the HF ground state is 
no longer a simple tensor product state even within the HF approximation and there is some nontrivial spin/valley 
entanglement. To facilitate subsequent discussion we summarize the results of Ref. [29] in Table [I] The analytical 
expressions for the mean-field levels were presented in Ref. E9] for the F/CAF cases. A straightforward calculation 
allows direct extension to the KD and the CDW configuration. Examples for the SP energy levels £±± obtained 
within this approximation for the different phases are plotted in Fig. [g] for various choices of the couplings. 

The ansatz that we introduced in Secs. [TT] and |IH] is able to describe spatial dependence of the order and also to 
capture spin/valley entanglement. 

We obtain the set of parameters {a(R), (3(R),x{R)} characterizing the GS | g), by minimizing the total energy E to t 
of Eq. (15). From these space-dependent parameters, we construct the corresponding density matrix p g = gg 1 via 


Eq. (|T 8 j) , which in turn allows us to reconstruct and diagonalize h HF of Eq. (|21[). The analysis is repeated for every 


point in space, thereby yielding the spatial behavior as a function of the distance from the edge. The results for the 
spectra £.; established near the edge for different phases in the bulk are shown in Figs. [T] and [ 8 } We have chosen the 
same parameters as in the plots of Figs. [3] and [5j For u z > 0 in Fig. [TJ we explore the upper half plane, whereas for 
u z < 0 (Fig. [8} , the evolution of the bulk phases of the lower half-plane is displayed. In all figures the thick colorful 
lines show the outcome of our numerical studies. The thin black lines represent the analytical results for the SP energy 
levels for the phase established in the bulk at the particular system parameters shown, as given in Table [TJ Figure [7] 
illustrates the evolution of the edge SP levels for a transverse coupling energy of u z = 5Ez and different values of the 
perpendicular coupling itj_, for which the bulk phase configuration passes from a CAF phase at u± = —Ez (green 
lines, upper left panel) to a F phase at u± = —0.2Ez,u± = 0.5Ez , and u_l = 1.5 Ez, respectively (blue spectra in the 
upper right and the two lower panels). In general, the spectra in all four panels show the following behavior : two flat 
energy levels, separated by the gap A £buiki are present in the bulk, both two-fold degenerate and they split into four 
branches when approaching the edge. The two intermediate levels, being labeled £2 and £ 3 , first bend towards each 
other, establishing the minimum energy gap A £ e d ge < A.£b u ik, before, even closer to the edge, the two lowest and the 
two highest levels £i ,£2 and £ 3 , £4 are driven apart in two parallel pairs, respectively. 

In Fig. [ 8 j we show the spectra corresponding to the phases of the lower half plane : at u z = —2, we display the 
energy levels at u_l = —0.6 Ez (red lines in the left panel), for which the bulk is CDW as well as for u± = 1.2 Ez 
(blue lines in the right panel), where the bulk is F. The behavior of the SP energy levels in a CDW bulk phase 
qualitatively differs from the situation of the CAF/F phase described above. In the left panel of Fig. [ 8 ] there are 
four non-degenerate levels in the bulk. In contrast to the levels of the CAF/F case, they do not bend towards each 
other and there is no minimum energy induced by the edge behavior. Hence, we find that the minimum edge gap is 
equal to the bulk gap : A £ e d ge = Ae bulk■ Sufficiently close to the edge, the levels again form two parallel pairs. The 
spectra for the KD phase in the bulk are not shown because, as mentioned in Sec. [TV) this state does not undergo 
any significant evolution when approaching the edge. The spectra do not differ from the analytical prediction for e±± 
shown for A A F 0 D in the lower right panel of Fig. pi 
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CAF/F phase: 


a A = a A 0z 1 0 cr z + A A zx a z 0 <Jx 

with a Aq z = — |(mo + tt z + 2 ui) cos# and a A zx = — |(mo + « z — 2ux) sin#, 

e±± = ±^/[-Bfcin(p) ± (Ez — A Aoz) ]" + ( A A zx y, 

As e dge — 2 | A zx |) A Ebulk ~ U 0 + Uz 2?iX, A £b u lk — 2|F/^ Aqz | ■ 


CDW/KD phase: 


a A = a A x0 ct x 01+ A A y oa y ®1+ a A z0 ct z 0 1 

with A X; yo — g (uot x ,y u z t X) y 4 Uj_t x ^ y ') and A z q — 2(^0 G 3 u z tz 2uj_G), 

which implies ^A^q 014 = — |(m 0 — 3u z — 2u±) and A A x0 D = — |(u 0 — u z — 4ux), 

4± W = ±Sz ± E kin (p) 2 + (- 4 A'?q DW ') 2 , and eg? = ±Fz ± [B w „(p) - a A“], 

Aebuik = 2\Ez — A A z @/ x p 11. 


TABLE I: Analytical expressions for the mean-field potential of the symmetry-breaking terms, A A, the eigenvalues of the full 
mean-field Hamiltonian, e±±, and the minimum gaps in the bulk and at the edge, A Sbuik/edge- We denote by s and t the SP 
spin and isospin configuration of the two electrons per orbital and 6 describes the canting angle between the two spins. In all 
the formulas for the mean-field potentials, we dropped a constant term — | (uo + 2mx +«.)101. These analytical results have 
been obtained within the approximation that the bulk phase does not change as a function of space when approaching the edge. 


When comparing this behavior to the analytical results as given in Table |T] (black lines), deep in the bulk we find 
that all curves coincide as they should. Furthermore, from the discussion in Sec. |IVB[ we know that there is no spin- 
valley entanglement in the bulk, i.e., the bulk states indeed are of separable product form. Yet, significant deviations 
between the numerical results capturing the full GS properties and the analytical curves from the simplified treatment 
are observed when moving closer to the edges, where the GS spin/isospin configuration starts to deviate from the 
bulk phase; cf. Fig. [3j The SP energies £j have kinks whenever the underlying spin/isospin texture changes and 
exhibit qualitatively different behavior in the different texture regimes. Thus, the emergence of different spin/isospin 
configurations due to the edge potential when approaching the edges directly translates into the SP spectra leading 
to a complex energy structure as a function of space. 


B. Single Particle Level Crossings in the different Texture Regimes 

The energy levels of the SP ground and excited states show a complex structure as a function of space depending 
on the spatial changes of the spin and isospin texture when approaching the graphene edge. In particular, in some 
configurations, the SP spectra exhibit a finite gap, whereas for other system parameters, the edge states are gapless 
when the SP levels cross. In this section, we investigate the crossings between SP energy levels which leads to gapless 
edge excitations. We first discuss the properties of the edge gap and its behavior when approaching the critical values 
where it closes. Then, we explain the number of allowed crossing points between energy levels and the connection to 
the symmetry properties of the underlying spin/isospin texture phases. The spatial variation of the order parameters 
has a direct impact on the overall shape of the dispersion of edge modes. This is readily seen in Figs. [7] and [ 8 j where 
we plot the dispersions from our calculation including edge effects and results from a similar calculation using only 
bulk values without spatial variation. In order to investigate the closure of the edge gap A e e d ge as a function of the 
ratio ux/F-Z) we evaluate the size of the minimum gap in the SP spectra for various system parameters. We choose 
the same values for the anisotropy energies as in Sec. |IV| by fixing u z = 5Ez,u z = 2 Ez and u z = —2 Ez and varying 
u± so that we access all bulk phases in the GS phase diagram. The resulting curves A £ e d g e{u±/Ez) are shown in 
Fig- [9} We find that the size of the edge gap Ae edge is a strictly monotonic decreasing function of u_\_lE z for all 
values of u z . When the bulk is KD or CDW the flat bulk SP levels split further apart when approaching the edge 
so that the minimum gap in the spectrum is equal to the bulk gap As edge = A Sbuik- In these two cases we find 
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FIG. 7: Spatial behavior of the energetic SP spectra in the presence of a boundary for positive transverse coupling u z = 5Ez 
and different values of the perpendicular coupling u± : CAF bulk phase at u± = —Ez in the upper left panel (green lines), 
F bulk phase at u± = —0.2 Ez, u± = 0.5 Ez, and u± = 1.5 Ez in the remaining panels, respectively (blue lines). Thick, 
colorful lines show our numerical results. Here, different line shapes distinguish between the different single-particle energy 
levels £i < £2 < £3 < £4. Thin, black lines compare to the analytical formulas for e±±(R) listed in Table [I] for the different 
phases, respectively, in which no modulation of the underlying spin/isospin texture towards the edge is taken into account. We 
see two bulk levels, being two-fold degenerate and separated by a gap of width Aebulk in the bulk, split into four branches when 
approaching the edge. The branches exhibit kinks and regimes of different behavior corresponding to the transitions between 
different spin and isospin phases during the evolution from the bulk to the edge (see Fig. [3b. The edge gap between the two 
intermediate levels £2 and £3 (dotted and dashed lines, respectively) Ae e d g e = min (|£3 — £ 2 |]remains finite over a certain range 
of ui, reducing gradually as u± increases until it finally closes completely. The lower right panel shows a configuration were 
the levels cross and form gapless edge states. Special attention should be paid to the upper right panel and the lower left panel 
in which the numerical results show configurations with a F phase in the bulk were finite edge gaps A £ e dge ^ 0 remain, whereas 
the analytical curves cross as soon as the bulk passes into an F phase. The behavior of the Ae e d g e as a function of u± is studied 
further in Fig. [9] 



ms. ... 1. r. ix... 1.... 1.... 1.... 1.... 1.... 1.... 1 

3 4 R 5 6 3 4 R 5 6 


FIG. 8 : Same as Fig. [7] for u z = —2Ez and different Mj_, hence displaying the SP spectra for the bulk phase regimes of the 
lower half plane of the GS phase diagram. Left panel: CDW bulk phase at u± = —0.6 Ez (red lines). Right panel: F bulk 
phase at u± = —1.3 Ez (blue lines). Thin, black lines compare to the analytical formulas given in Table [I] for the different 
phases, respectively, in which the evolution of the bulk’s phase towards the edge is neglected. For the CDW wave phase in the 
left panel we observe four non-degenerate SP levels in the bulk; furthermore, as these levels do not bend towards each other 
but disperse when approaching the edge, the minimum energy gap to SP excitations is given by the bulk gap Ae e d g e = A Ebuik- 


that the bulk gap is a linear function of the perpendicular coupling energy : A.£ 1 ^j ) l / CDW oc u±/Ez . The numerical 
results in Fig. |9]follow exactly the analytical prediction given in Table]!]: At u z = —2 Ez, we find A = —4u_l and 
A ^buik' = — 2u_l +4 Ez, whereas the KD edge gap at u z = 2 Ez behaves as Ae^jfc = — 4u_l — 4 Ez- These analytical 
curves are plotted in Fig. [9] as dotted lines for comparison (they are shifted by a constant offset with respect to the 
numerical results for better visibility). As a consequence of this linear behavior in itj_, for couplings favoring KD 
or CDW in the bulk, at the system parameters chosen in Fig. [9] there is always a non-zero gap to SP excitations. 
When the GS in the bulk is a CAF or a F phase, the SP spectra now bend towards each other and therefore exhibit 
a minimum energy gap A £ e d ge near the edge which is smaller than the bulk gap. 
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FIG. 9: Behavior of the gap A £ edge in the SP spectra close to the edge. We use three couplings : u z = —2 Ez, u z = 2 E z , and 
u z = 5 Ez (filled, shaded and empty symbols, respectively). Different colors and symbols stand for different bulk phases: KD 
(grey triangles), CDW (red diamonds), CAF (green squares), and F phase (blue circles). Dashed, colored lines connect the 
data points as a guide to the eye. The dotted, black lines represent the behavior of the data in the linear regimes (they are 
shifted by a constant offset with respect to the curves for better visibility). Gray vertical lines indicate the critical values for 
bulk phase transitions. For all phases the gap A £ edge monotonically decreases as u± grows until it finally closes in the regime 
where the bulk is in an F phase. The inset is a close-up on how the blue lines smoothly approach Ae = 0. The transitions from 
A £ e dge 7 ^ 0 to AE edge = 0 take place at u± « 0.3 E z , u± « E z , u± « 1.625Uz. 


For u _l < —Ez /2 where the bulk is a CAF phase (green squares), the spectrum always exhibits an non-zero edge 
gap which is almost linear as a function of the perpendicular coupling. At values u± > —Ez/ 2, hence for a F bulk, the 
shape of the spectrum changes qualitatively and the bulk gap closes in a non-linear way, asymptotically approaching 
zero at sufficiently large values of uj_. For the transverse couplings chosen in Fig. [9j u z = 5Ez, u z = 2 Ez, and 
u z = —2 E z , the edge gap Ae edge closes at u± ~ E z , u± « 0.3E Z , and u± « 1.6E Z , respectively. All these 
values leads to a bulk which is ferromagnetic. So the prediction for the gap closure point clearly differs from the 
value u± = which can be read from the bulk phase diagram in Fig. [ 2 ] This is due to the changes of the 

spin/isospin configuration of the GS induced by the effective edge potential as we approach the boundary; cf. Fig. [ 3 ] 
of Sec. |IV| Indeed the system does not remain in a F phase configuration all the way from the bulk to the edge. 
During its transition into a KD phase close to the boundary, there is an intermediate regime with non-zero spin-valley 
entanglement and simultaneous canting of both spin and isospin. Hence, in this transition regime there is no a priori 
justification for of Table i] to yield a correct description of the edge gap. 

From this analysis of the gaps 01 the SP spectra we can draw the following conclusion : when the bulk is CDW, KD 
or CAF phase, the SP energy levels always have non-zero gaps. However for a bulk F phase one may have gapped 
or gapless spectra. We note that ignoring the spatial variation of the trial HF state leads to qualitatively different 
result^. 


C. Number of level crossings 

We now discuss in more detail the number of crossings of the HF single-particle states. In the F phase with 
dispersion e^ ± in Table |lj the intersecting levels and e F /_ cross exactly once as their slope is given by the slope 
of the kinetic energy term and they are monotonic functions of the spatial coordinate. For several system parameters, 
such as in the spectrum in the lower left panel of Fig. [7] at u z = 5Ez,u± = 1 .5Ez as well as in close ups shown in 
Fig. [lOj we observe multiple crossings. We first discuss the occurrence of multiple crossings and the relation with 
the underlying spin/valley texture. The number of crossings is governed by the symmetries of the HF Hamiltonian 
and the magnitude of the HF self-consistent potentials. After discussing the different phases separately, we apply the 
insights to the edge state structure described in Sec. |IV[ where the GS phase changes as a function of space when 
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FIG. 10: Close-up on the SP spectra in the regime with multiple crossings for u z = 2 Ez- Left side: u± = 1.2 E z , right side: 
u± = 3Ez . The blue lines show the SP energy levels £,. The different background colors mark the spatial regimes in which the 
GS establishes different spin/isospin textures: there is a F phase (blue region). In the intermediate (white) region, the system 
undergoes a transition before it finally ends up in a KD phase (yellow region). We observe the crossings of the SP levels to 
occur in regions of different spin/isospin phases - their origins lie in the different symmetries of the F and KD phases. 


approaching the edge from the bulk. We first discuss the case of the CAF/F transition. We rewrite the mean-field 
Hamiltonian of Eq. (21) involving the CAF/F mean-field potential given in Table [T] by decomposing it into four 2x2 
matrices as : 


HF 


0 ) = 


7i 72 
73 74 


(24) 


where the respective entries for the CAF/F phase are given by : 

71 = A A zx a x - ( E z - j 4 A 0z ) ct,, 

72 = 73 = -Ekinip) 1 , 

73 = - ' 4 A zx cr x - {Ez - A A 0z ) a z , 

with 4 A S2 , and a Aq z defined for the CAF/F phase in Table hi The size of the gap is therefore governed by the first 


(25) 


off-diagonal coupling matrix elements A A 


A A CAF/F 


ZX 

HF 


If a A < ^ AF ' F 0, as is the case for any non-zero canting angle 
0/0, the eigenvalues of the Hamiltonian h caf/f exhibit the characteristic behavior of avoided crossings. The SP 

levels are allowed to cross only for A A z ' 4F/ ” F = 0 at 9 = 0, i.e., in the F phase. In the bulk, i.e., at Eki n = 0, all 
values of the coupling energies u z and u± allowed for the F phase yield the same ordering of the SP energy levels 
e± 4 l = e±±{Ekm = 0), independently of the sign or the modulus of Aq AF ^ F : e F ’® = e F ’° < 0 < s F '°_ = e F +- Hence, 
there is only one possible scenario of level crossings when approaching the boundary as the increasing edge potential 
is driving the SP levels away from their bulk values. This leads to exactly one crossing of the levels e+_ and £ F _, 
shown by the blue, dashed lines in the upper left panel of Fig. [bj We can perform the same analysis for CDW or 
KD phases. Again, we rewrite the corresponding HF Hamiltonians of Eq. (21) with the potentials A A CF ^J KD from 
Table [I] and we find the respective entries for the CDW : 


whereas for the KD phase we find : 


71 = —Ez cr z + a A z0 1, 

72 = 73 = - Ekinip) t, 
74 = -E z cr z - A A, 0 1, 


7 1 = 7 4 = E z a z , 

72 = 73 = ( AA xO D - Ekinip )) 1 - 


(26) 


( 27 ) 
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The Hamiltonians for the CDW phase and the KD phase thus turn out to have higher symmetry than in the CAF 
phase: In h C£W and h KD , all entries of the two first off-diagonals as well as of the antidiagonal are zero. Pairwise 
degeneracy of the corresponding eigenvalues, i.e., crossings between the SP energy levels, is now allowed. Note that, 
unlike the transition from a CAF to a F phase, all other transitions do not correspond to smooth transitions. In these 
cases, a transition between phases goes along with an abrupt change of the symmetry properties of the spin/isospin 
configuration of the GS and the corresponding Hamiltonian. 

We now discuss the different possible scenarios of SP level crossings in CDW and KD phases. The SP energy 
levels of the CDW s^ w in Table 0 are independent of the sign of A A^ 0 DA . Different orderings of the bulk levels 
at Ekin = 0 may, however, appear depending on the modulus of A A?,f w : for | a A^ 0 dw '| > Ez, the 
bulk states are ordered as ,0 < In this case, when approaching the boundary, 

the kinetic energy potential drives the positive and the negative energy states further apart from each other such 
that they do not cross. In the case where | A A^ ( f >n \ < Ez, however, the bulk states rather follow the hierarchy 
£ cnwfi < £ cdu, o < q < £ Cdw,o < £ Cdw,o ^is casej turning on the effective edge potential drives the levels 

£ CZJU ,o an j g’t-° towards each other and they cross at zero energy. These two different scenarios are depicted in 
the upper right panel of Fig. [bj where the red, solid lines show the levels e±± w from Table[l]at A A ^ 0 DW = 2 Ez > Ez 
and the black, dashed lines show the spectrum for A A^ W = 0.3Ez < Ez- The latter case, ' 4 A(? 0 DM ' \ < Ez, 
however, is prohibited by the conditions imposed on the couplings u z and zi_L in order for the system to establish 
a CDW phase in the bulk. Requiring u z < u± and u z < —Ez — u± will always force | A A^ 0 nw | > Ez- Therefore, 
treating the system as having a stable CDW phase in the bulk and all the way to the edge will never lead to any 
crossings of the SP edge levels. Turning to the more important case of the KD phase, the situation becomes even 
richer. Here, depending on the sign and the modulus of A A^ 0 D , four different SP level orderings in the bulk and four 
resulting crossing scenarios may appear. For I^A^I < Ez, if A A™ > 0, there is one level crossing at zero energy 
and two additional crossings above and below the zero energy line, respectively, whereas for negative a A^ 0 d , only 
one crossing at zero energy is present. The case | A A^j°| > Ez can lead to four crossings, two at zero energy plus one 
above and one below, respectively, if A ! ^ 0 D > 0 , whereas for A A^ 0 D < 0 , the four levels do not cross. 

Examples of the four different cases are shown in the lower panels of Fig. | 6 j where the lower left panel shows 
the possible situations for A A^ 0 n \ < Ez (solid, orange lines for A A£ 0 d = —0 -3Ez < 0 and black dashed lines for 
A A« D = 0.7Ez > 0), whereas the right panel displays the corresponding spectra for | a A£ 0 d \ > Ez (here, the solid, 
orange lines are for A A A 0 D = — 1.5Ez < 0 and black dashed lines for A A^ 0 D = 2 E z > 0). Note that, just as in the 
case of the CDW phase, not all these cases are allowed by the restrictions on the parameter range for a KD bulk : 
requiring the couplings u z and u± to fulfill the relations u_ j_ < u z and u z < ^— u± always implies a A^ 0 d > Ez- 
Therefore, again, all cases including possible crossings between edge levels are ruled out for a system with a KD phase 
in the bulk. 


This simple picture drawn for constant order parameters changes when considering the electronic GS structure 
described in Sec. IV Indeed the GS spin/isospin texture deviates from the bulk phase when moving towards the 
edge as a consequence of the growing edge potential. Close enough to the edge the system is always driven into a 
KD phase. Hence when moving sufficiently close to the edge the GS does becomes KD-ordered, even though the 
system parameters u z and u± do not allow KD order in the bulk. Two examples are shown in the close-ups in 
Fig. [lO] Parameters in both panels are chosen such that the bulk system at E^ in = 0 is in a F phase. When moving 
towards the edge the energy levels evolve according to e± ± of Table [I] (corresponding to the evolution within the blue 
region). A first crossing between the intermediate levels occurs, as predicted by the analysis of the F phase energy 
levels. After the transition region (left white), the GS becomes KD (marked by the yellow shading). However, the 
system parameters do not force A A^ 0 D > Ez : in the left panel of Fig. 10 we find A A 
right panel we have A A £ 0 n = —7 Ez- Therefore the energy levels now evolve according to e 


_KD 
z ±± 


of Table [T] ir 


the 

case A A^ 0 D < 0, | a A£ 0 d \ < Ez- As a consequence, in this regime, one more level crossing may occur. Hence, the 
appearance of several crossings of the SP energy levels in the numerical spectra as in the lower right panel of Fig. [7] 
and in Fig. [TO] can be explained combining the insight of Sec. |IV| that any bulk phase by the edge potential always is 
driven into a KD phase close to the boundary, with the understanding of the possible behavior of depending on 
the value of A A^ 0 D as a function o the coupling energies u z and u_l. The SP energy levels describing the numerical 
results of Figs. mu and [ 10 ] can be summarized as : 


r 4“± fe (-R) 

£±±(T?) = ( unknown 


for R> R 2 , 
for Ri < R < R 2 , 
for R < Ri, 


(28) 


where e±±(R) denotes the level spectra for the bulk phase established at a given choice of system parameters and R 2 
and R\ label the inner and outer limits in space of the domain wall, for which there is no simple analytic expression. 
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FIG. 11: Evolution of the single-electron spin and isospin components s z (i) and t x (i) of the SP eigenstates from the bulk 
towards the edge for the two anisotropy energies m_l = — Ez (CAF bulk GS , left panels with green lines) and ui = —0.2 Ez (F 
bulk GS, right panels with blue lines). Different line shapes distinguish between the four SP energy levels £i < £2 < £3 < £ 4 . 
Green/blue lines correspond to the observables of the two lowest-lying states which are occupied orbitals in the HF GS. Gray 
lines indicate the behavior of the higher-lying SP states. Arrows show the behavior of the spin and isospin polarizations. The 
second and third levels \2) and |3) are oppositely polarized in spin and isospin at the edge. In all plots we set u z = 5Ez- 


The evolution of ( R ) is no longer limited to the non-crossing behavior imposed for a bulk KD phase but it can 

exhibit any of the shapes drawn in the two lower panels of Fig. [6] Which of these curves describes the KD-like 
evolution of the edge states correctly is determined by the system parameters and u z that govern the bulk texture 
phase. From the analysis of the number of SP level crossings we hence learn that, in principle, by choosing appropriate 
values of u z and itj_, SP energy levels can have zero, one, two, or even three crossings at zero energy. Among these 
crossings, one is due to the symmetry properties of the bulk F phase. The remaining crossings appear in the KD 
phase close to the boundary, which in this regime shows novel properties not present for a KD phase in the bulk. We 
note that these crossings occur at different distances from the edge - at the distance where the corresponding KD 
phase SP levels for a certain A A^ 0 D cross, it is necessary for the system already to have evolved from the bulk phase 
into the KD edge phase in order for the additional SP level crossings to occur. This is the reason why we do not 
see any crossings in the SP spectrum shown in the lower left panel of Fig. [8] where the bulk is in a CDW phase : at 
the distance R cr oss where the KD-like levels near the edge would cross, the system still behaves according to its bulk 
CDW phase. In this case, the crossing is thus prevented by the fact that the crossing point lies outside the KD region : 
Rcross > R 2 • Nevertheless, a situation in which the bulk is a CDW but the SP edge states are gapless due to crossings 
of the KD-like levels close to the boundary is not forbidden by the underlying symmetry principles as the restrictions 
for the coupling energies of the CDW bulk phase allow negative values of A A ^ {) n . The exact distances from the edge 
or Rcross which define the points of crossing, involve the explicit form of the kinetic energy Eki n (R) as they 
are determined by the eventual dominance of the kinetic energy. Numerical values for i?i,I? 2 i or R cr oss therefore 
strongly depend on the model potential chosen for Ek in (R). This is not true, however, for the answer to the question 
of whether crossings are allowed or not since the values of ^ A are determined generically by the system parameters 
u± and u z . 


D. Properties of the underlying SP States 

In order to obtain a better understanding of the nature of the excited states we analyze the properties of the single¬ 
electron states. We compute the spin and isospin components s z (i) = \{i\a z \i) and t x (i) = \{i\r x \i) and display the 
results as a function of space in Fig. El 

The evolution of the observables can be summarized as follows : 
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where arrows represent schematically the spin and isospin vectors. 

The HF GS |G) is built from Slater determinants of |1) and |2) as the eigenstates of the two lowest lying branches 
£1 and £2 in Fig. [7] The lowest energy excitations correspond to single-particle excitation from the second to the third 
level £2 —> £ 3 . These two states have oppositely polarized spin and isospin components. The closing of the gap A £ e d ge 
between the second and the third SP level in Fig. |TT] hence is a transition from insulating to conducting behavior 
with the counterpropagating current-carrying edge states exhibiting opposite spin and isospin polarizations. This is 
the behavior of gapless helical edge states of a QSH stated. 


VI. CONCLUSION 

The SU(4) QH ferromagnetism leads to a highly degenerate manifold of ground states for neutral graphene. This 
degeneracy is lifted by small lattice-scale anisotropies and there is a competition between phases with several types of 
order. This competition is affected notably by the substrate supporting the graphene sample. We have used a simple 
model of these anisotropies to study the edge properties of neutral graphene by means of a HF approach. The sharp 
atomic edge is then described by an effective field in valley space which modifies the competition between phases. 
Ultimately, whatever the bulk order, the system is in a KD phase close enough to the edge. In the transition region 
between this edge order and the bulk order, we have obtained evidence for an intermediate regime with spin/valley 
entanglement. In this regime there is a nontrivial change of the single-particle spectrum. We find that the number 
of levels crossing the Fermi energy can be varied by changing the parameter u±/Ez . This means that there are 
metal-insulator transitions when tilting the magnetic field. This is consistent with the experimental findings of Young 
et all 22 . If we adopt the estimates for the approximate magnetic field dependencies of Ez and u_l stated in Refs. (7] 
and 1251 as E Z {B) k, 0.7B [T]K and u±(B±) ss 1 — 10Rj_[T]K, where B denotes the total magnetic field and B± its 
component perpendicular to the device plane, the values for the parameters stated in Ref. [2H suggest that the authors 
were able to experimentally tune the ratio u±/Ez roughly in a range from -13 to -0.5. The picture we obtain is more 
complex than that obtained by assuming that the order does not persist up to the edge®( Notably the occurrence of 
the metal-insulator transition, while it sets constraints on the microscopic parameters, does not imply that the bulk 
is CAF ordered. The observation of a conductance G « 2 e 2 /h which corresponds to two conducting channels, i.e., 
to one single level crossing, has two possible explanations : either the bulk is in a F phase leading to one crossing 
unaffected by the KD edge regime, or the bulk has noncrossing SP levels, but the crossing occurs in the KD regime 
close to the edge. Furthermore, our results suggest that the observation of exactly one crossing only corresponds to 
a limited parameter range. Varying the anisotropy parameters may lead to the observation of conductance values of 
higher multiples of two, corresponding to several crossings in the SP edge spectrum. 

Of course there are obvious limitations of our theoretical approach : we apply a perturbative treatment of the edge, 
whose validity is limited to a certain range [cf. the discussion following Eq. <[s|)]. The appearance of a KD phase in 
the vicinity of the edge is a direct consequence of treating the effective edge potential perturbatively. Furthermore, in 
our calculations we assume the anisotropy energies u± and u z to remain constant at their bulk values [as we discuss 
after introducing Eq. 0 ]- This approximation certainly becomes less exact as we approach the boundary. Also we 
have neglected the exchange energy effects that will create textures in the charge-carrying states. 

In conclusion, in this paper we have studied the influence of an edge on the v = 0 QH state in monolayer graphene. 
We found that the effective edge potential induces a change of the GS spin/isospin texture. During this evolution, novel 
phases are observed, involving simultaneous canting of spin and isospin as well as non-zero spin/valley entanglement. 
Phases of this kind are not present in the bulk. Furthermore, we analyzed how this spatial evolution changes the SP 
excited states. Here we have shown that, as a consequence of the spatial modulation of the underlying spin/isospin 
texture, the direct correspondence between the conductance properties of the edge states and the bulk phase is lost. 
The transport properties are governed by either zero, one, or multiple SP level crossings. The analysis of the SP 
eigenstates shows that the lowest SP excitation describes counterpropagating helical edge states carrying opposite 
spin and isospin polarizations. 
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